#### carregando pacotes e scripts ####
#====================================#
knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE
)
options(outDec=",",big.mark=".")
set.seed(3)
library(flextable)
library(tidyverse)
library(randomForest)
library(naivebayes)
library(plotly)
library(MLmetrics)
set_flextable_defaults(scroll=list())
source("src/functions.R")
source("src/texts.R")
source("src/setup.R")As variáveis mais indicadas para alimentar um modelo de previsão já foram discutidas anteriormente. Agora estas variáveis serão agregadas por internação e separadas em dois grupos: um para treinamento dos modelos e outro para avaliar a qualidade das previsões. O grupo de treinamento será uma amostra aleatória contendo 80% dos dados disponíveis, enquanto que os dados de avaliação serão os demais 20% dos dados que não foram incluídos no grupo de treinamento.
Os principais negócios interessados no valor deste conjunto de dados seriam as empresas seguradoras de saúde, estas estariam principalmente em prever se haveria necessidade de transferir algum recurso para a fornecedora de saúde (Hospitais e Clínicas, por exemplo).
Podemos melhorar a qualidade de vida de tais empresas ajudando a reduzir as incertezas sobre o momento em que serão acionadas. Para isso, vamos adicionar uma variável-alvo que represente o evento de acionamento do seguro.
# Criando variavel target
agg_df$sinistro <- agg_df$VL_ITEM_PAGO_FORNECEDOR > 0
# Obtendo informações sobre o balanceamento dos dados
QTD_ACIONAMENTOS <- sum(agg_df$sinistro)
QTD_CASOS <- nrow(agg_df)A próxima etapa é escolher as variáveis que podem ajudar nosso modelo preditivo a chegar na resposta esperada com tempestividade. Mas apenas depois que obtivermos um conjunto de dados balanceados, pois apenas 3.84% das ocorreências registradas apresentam acionamento de seguro, para evitar vieses, o efeito da simples quantidade de observações não pode ser maior que a presença das variáveis explicativas.
# Obtendo indices de `agg_df`
df_index <- row.names(agg_df) %>% as.numeric()
# Obtendo indices das amostras balanceadas
### maximizando a quantidade de observações usando todos os casos positivos
negative_sample_index <- sample(
x=df_index[!agg_df$sinistro], size=QTD_ACIONAMENTOS)
positive_index <- df_index[agg_df$sinistro]
mdl_index <- c(negative_sample_index, positive_index)
# Selecionando variáveis úteis
mdl_cols <- c(
"NM_MODALIDADE", "UF_PRESTADOR", "FAIXA_ETARIA",
"ANO_MES_EVENTO", "sinistro")
# Criando novo dataframe com valores de "sinistro" balanceados
mdl_df <- agg_df[mdl_index, mdl_cols]Já aprendemos anteriormente como estas variáveis categóricas se relacionam com o acionamento do seguro, o objetivo desta engenharia de variáveis é simplificar a interpretação do modelo, mas principalmente reduzir a quantidade de features aumentando sua importância.
Começando pelo NM_MODALIDADE, sabemos que
“Cooperativa médica” e “Medicina de Grupo” são os tipos de prestadoras
com uma maior proporção de acionamentos de seguros, portanto vamos
reduzir esta variável para apenas dois níveis, um para o caso de
participar de alguma delas, outro para o caso contrário.
temp <- mdl_df$NM_MODALIDADE
temp[temp %in% c("Cooperativa Médica", "Medicina De Grupo")] <- "Cooperativa"
temp[temp %in% c("Seguradora Especializada Em Saúde",
"Filantropia", "Autogestão", NA)] <- "Não-cooperativa"
mdl_df["modalidade"] <- factor(temp)Percebemos também que em UF_PRESTADOR, os
estados do nordeste se tornam mais presentes em relação aos demais,
portanto reduzir os estados às suas regiões já deve ser suficiente.
temp <- mdl_df$UF_PRESTADOR
temp[temp %in% c("MA", "PI", "CE", "RN", "PB", "PE", "AL", "SE", "BA")] <- "Nordeste"
temp[temp %in% c("AC", "AP", "AM", "PA", "RO", "RR", "TO")] <- "Norte"
temp[temp %in% c("GO", "MT", "MS", "DF")] <- "Centro-Oeste"
temp[temp %in% c("ES", "RJ", "SP", "MG")] <- "Sudeste"
temp[temp %in% c("PR", "RS", "SC")] <- "Sul"
mdl_df["regiao"] <- factor(temp)Embora os mais idosos estejam entre os mais atendidos pelos
fornecedores, os adultos estão entre os maiores acionadores de seguro,
vamos reduzir esta quantidade excessiva de
FAIXA_ETARIA em apenas 3 (crianças, jovens/adultos
e idosos), os casos não identificados formarão um quarto grupo.
temp <- mdl_df$FAIXA_ETARIA
temp[temp %in% c("Não identificado", NA)] <- "Não Identificado"
temp[temp %in% c("<1", "1 a 4", "5 a 9")] <- "Crianças"
temp[temp %in% c("10 a 14", "15 a 19", "20 a 29",
"30 a 39", "40 a 49", "50 a 59")] <- "Jovens/Adultos"
temp[temp %in% c("60 a 69", "70 a 79", "80 ou mais")] <- "Idosos"
mdl_df["idade"] <- factor(temp)Por fim, percebemos também que os seguros são acionados menos vezes
nos primeiros três meses do ano, em detrimento dos demais. Vamos agregar
os dados em ANO_MES_EVENTO por trimestres,
independentemente do ano.
temp <- mdl_df$ANO_MES_EVENTO |> format("%m") |> as.integer()
temp[temp %in% 1:3] <- "Primeiro"
temp[temp %in% 4:6] <- "Segundo"
temp[temp %in% 7:9] <- "Terceiro"
temp[temp %in% 10:12] <- "Quarto"
mdl_df["trimestre"] <- factor(temp)Vamos treinar alguns modelos de regressão logística para obter a melhor combinação de variáveis no trabalho de identificar a se um seguro será ou não acionado.
logit_model <- function(formula, data)
glm(formula=formula, data=data, family=binomial(link="logit"))# exemplos de modelos testados
mdl_v1 <- logit_model(sinistro ~ modalidade*idade -1, data=mdl_df)
mdl_v2 <- logit_model(sinistro ~ modalidade * trimestre, data=mdl_df)
mdl_v3 <- logit_model(sinistro ~ idade * trimestre * regiao, data=mdl_df)
mdl_v4 <- logit_model(sinistro ~ modalidade + idade * regiao, data=mdl_df)
mdl_v5 <- logit_model(sinistro ~ idade * modalidade * regiao, data=mdl_df)
mdl_v6 <- logit_model(sinistro ~ idade + modalidade * regiao, data=mdl_df)
mdl_v7 <- logit_model(sinistro ~ trimestre + modalidade * regiao, data=mdl_df)
mdl_v8 <- logit_model(sinistro ~ trimestre + idade + modalidade * regiao, data=mdl_df)
mdl_v9 <- logit_model(sinistro ~ idade + trimestre * modalidade * regiao, data=mdl_df)Depois de passar um tempo testando as variáveis e suas interações, as
melhores combinações são aquelas em que o trimestre e a
idade não interagem com as demais variáveis, entre as
modelagens listadas acima, as de número 6, 7 e 8 mostraram coeficientes
estatisticamente mais significantes que os demais, por obedecerem esta
condição.
Usando o Critério de Informação de Akaike (AIC) como desempate, a
formulação do mdl_v8 se mostrou a mais eficiente. Vamos
considerar esta modelagem nas análises adiante, mas antes, vamos ver
algumas estatísticas em relação à este modelo.
form1 <- sinistro ~ trimestre + idade + modalidade * regiao
mdl_lg <- logit_model(form1, data=mdl_df)
summary(mdl_lg)##
## Call:
## glm(formula = formula, family = binomial(link = "logit"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.87573 0.17577 -4.982 6.28e-07
## trimestreQuarto 0.57123 0.08622 6.625 3.46e-11
## trimestreSegundo 0.35421 0.08661 4.090 4.32e-05
## trimestreTerceiro 0.45544 0.08612 5.288 1.24e-07
## idadeIdosos 0.71699 0.13319 5.383 7.32e-08
## idadeJovens/Adultos 0.84409 0.12808 6.591 4.38e-11
## idadeNão Identificado 0.90338 0.22238 4.062 4.86e-05
## modalidadeNão-cooperativa -2.30599 0.29274 -7.877 3.35e-15
## regiaoNordeste 1.89774 0.13970 13.584 < 2e-16
## regiaoNorte 1.74420 0.20055 8.697 < 2e-16
## regiaoSudeste -0.81231 0.12773 -6.360 2.02e-10
## regiaoSul -0.55502 0.14082 -3.941 8.10e-05
## modalidadeNão-cooperativa:regiaoNordeste -2.08671 0.36754 -5.678 1.37e-08
## modalidadeNão-cooperativa:regiaoNorte -15.16443 205.34828 -0.074 0.9411
## modalidadeNão-cooperativa:regiaoSudeste 0.58922 0.31486 1.871 0.0613
## modalidadeNão-cooperativa:regiaoSul 1.38438 0.35080 3.946 7.94e-05
##
## (Intercept) ***
## trimestreQuarto ***
## trimestreSegundo ***
## trimestreTerceiro ***
## idadeIdosos ***
## idadeJovens/Adultos ***
## idadeNão Identificado ***
## modalidadeNão-cooperativa ***
## regiaoNordeste ***
## regiaoNorte ***
## regiaoSudeste ***
## regiaoSul ***
## modalidadeNão-cooperativa:regiaoNordeste ***
## modalidadeNão-cooperativa:regiaoNorte
## modalidadeNão-cooperativa:regiaoSudeste .
## modalidadeNão-cooperativa:regiaoSul ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10485.9 on 7563 degrees of freedom
## Residual deviance: 7191.3 on 7548 degrees of freedom
## AIC: 7223.3
##
## Number of Fisher Scoring iterations: 14Vamos adicionar mais alguns modelos com a mesma configuração de variáveis, para este caso, vamos incluir os modelos “Naive Bayes” e “Random Forest”. Este último é conhecido na literatura por ser usado em modelagens de Churn e desempenhar bem, comparado com os demais que estamos usando.
form2 <- factor(sinistro) ~ trimestre + idade + modalidade * regiao
mdl_nb <- naive_bayes(form1, data=mdl_df)
mdl_rf <- randomForest(form2, data=mdl_df)A métrica de Akaike foi muito útil para escolher entre as regressões logísticas, mas não vai ser útil para comparar a regressão com os demais modelos, para isto, vamos computar a Área Sob a Curva ROC (AUC), esta curva delimita a relação de custo/recompensa de verdadeiro-positivo e falso-positivo para cada limiar de probabilidade de acionamento do seguro escolhida.
Para aumentar a quantidade de verdadeiro-positivos, um limiar de
probabilidade mais alto deve ser escolhido, mas ao custo de aumentar o
número de falso-positivos. Esta “taxa de troca” não é constante, por
isso se trata de uma curva, e quanto maior a àrea sob esta curva, melhor
o modelo se ajusta para toda combinação de limiares possíveis. Vamos
obter a curva ROC considerando dados em df, que representam
melhor as condições dos dados na vida real, para isto, vamos replicar o
feature engineering que fizemos antes:
df$sinistro <- df$VL_ITEM_PAGO_FORNECEDOR > 0
temp <- df$NM_MODALIDADE
temp[temp %in% c("Cooperativa Médica", "Medicina De Grupo")] <- "Cooperativa"
temp[temp %in% c("Seguradora Especializada Em Saúde",
"Filantropia", "Autogestão", NA)] <- "Não-cooperativa"
df["modalidade"] <- factor(temp)
temp <- df$UF_PRESTADOR
temp[temp %in% c("MA", "PI", "CE", "RN", "PB", "PE", "AL", "SE", "BA")] <- "Nordeste"
temp[temp %in% c("AC", "AP", "AM", "PA", "RO", "RR", "TO")] <- "Norte"
temp[temp %in% c("GO", "MT", "MS", "DF")] <- "Centro-Oeste"
temp[temp %in% c("ES", "RJ", "SP", "MG")] <- "Sudeste"
temp[temp %in% c("PR", "RS", "SC")] <- "Sul"
df["regiao"] <- factor(temp)
temp <- df$FAIXA_ETARIA
temp[temp %in% c("Não identificado", NA)] <- "Não Identificado"
temp[temp %in% c("<1", "1 a 4", "5 a 9")] <- "Crianças"
temp[temp %in% c("10 a 14", "15 a 19", "20 a 29",
"30 a 39", "40 a 49", "50 a 59")] <- "Jovens/Adultos"
temp[temp %in% c("60 a 69", "70 a 79", "80 ou mais")] <- "Idosos"
df["idade"] <- factor(temp)
temp <- df$ANO_MES_EVENTO |> as.Date() |> format("%m") |> as.integer()
temp[temp %in% 1:3] <- "Primeiro"
temp[temp %in% 4:6] <- "Segundo"
temp[temp %in% 7:9] <- "Terceiro"
temp[temp %in% 10:12] <- "Quarto"
df["trimestre"] <- factor(temp)Depois vamos guardar uma amostra destes dados para avaliar os modelos
somente depois de realizar todos os experimentos. assim teremos um
conjunto de dados que sabemos que nunca foi usado para treinamento de
modelos, para isso vamos usar uma amostra de 10% dos dados em
df:
var_cols <- c("sinistro", "modalidade", "regiao", "idade", "trimestre")
df_index <- 1:nrow(df)
adf_index <- sample(x=df_index, size=129451)
# dados de avaliação
adf <- df[adf_index, var_cols]
# dados de treinamento
tdf <- df[!(df_index %in% adf_index), var_cols]Agora podemos obter os valores das previsões de cada modelo, e colocá-las lado a lado em um data frame:
pred_df <- data.frame(
y=tdf$sinistro,
Regressao.Logistica=predict(mdl_lg, tdf),
Naive.Bayes=predict(mdl_nb, tdf, type="prob")[, "TRUE"],
Random.Forest=predict(mdl_rf, tdf, type="prob")[, "TRUE"]
)Vamos usar uma validação cruzada, seguindo o método
de Monte Carlo, mas preservando um conjunto balanceado de treino. O
código fonte deste algorítimo e mais informações sobre seu funcionamento
podem ser encontrados em functions.r.
Manter o equilíbrio da quantidade de observações de cada classe no treino gera dois problemas: 1- Diminui a quantidade de observações usadas no treino e 2- Desequilibra ainda mais os dados no teste. Mas estes custos não superam o benefício da remoção de viés no momento do treinamento, a maioria dos modelos pode simplesmente aprender que um evento ocorre mais que o outro, e desconsiderar o impacto das variáveis explicativas.
O método de Monte Carlo foi escolhido no lugar do K-Fold, pois este apresenta menor variância para as métricas, o que será um bom equiliíbrio, uma vez que já estamos controlando os modelos para viés. Assim estamos buscando um equilíbrio no trade-off conhecido entre viés e variância. A quantidade de 16 repetições foi escolhida devido a limitações computacionais.
cvdata <- mccv_data(
formulas=c(form1, form1, form2),
models=c(LogReg=logit_model, NBayes=naive_bayes, RndFst=randomForest),
types=c("response", "prob", "prob"),
dataset=tdf,
balanced_for="sinistro",
n_repeats=16)Agora que temos os dados de erros dos modelos, podemos obter métricas customizadas, não precisamos nos limitar ao que os pacotes do R oferecem por padrão. Para o propósito deste desafio, vamos usar uma métrica customizada que seja mais adequada possível para o problema que estamos enfrentando.
Como estamos usando um modelo em que as previsões são desbalanceadas, e em que o falso negativo (erro tipo 2) oferece mais custos de negócio que os falsos positivos (erro tipo 1), nossas métricas devem priorizar a capacidade do modelo de prever corretamente os casos positivos enquanto minimiza os erros tipo 2. Por este motivo, nossa métrica \(M\) incorpora o Valor Preditivo Positivo (ou Precision) e a Taxa de Falso Negativo:
\[ M_L = VPP_L - TFN_L \]
O valor previsto pelos modelos no caso de uma classificação binária é interpretado como a probablidade de ocorrência do caso postivo (acionamento do sinistro), há um limiar de 50% definido por padrão na maioria dos casos de previsão, onde a probabilidade maior que este valor é entendida como uma previsão de que o caso positivo acontecerá. Os modelos de classificação podem ser ajustados para capturar mais os casos positivos ou os casos negativos quando alteramos o limiar da probablidade prevista.
Esperamos que nosso modelo seja capaz de maximizar a proporção de valores positivos, para cada limiar \(L\), nós vamos obter o valor de \(M_L\), que pode oscilar entre 1 e -1, e quanto maior o valor obtido, melhor é o desempenho do modelo nas restrições estabelecidas.
Podemos obter um valor de “Area sob a curva” como se obtém com uma métrica conhecida chamada AUC. Mas como neste caso, nossa curva pode aparecer no lado negativo, é possível obter valores negativos desta métrica, podemos considerar um bom modelo aquele que apresenta valor positivo em sua “Área sob a curva”, sendo o modelo perfeito aquele com área igual à 1, este teria um bom desempenho .
for (i in seq_along(cvdata)) {
mdl <- names(cvdata[i])[1]
if (i == 1)
mauc <- MAUC(cvdata[[i]]$y_actual, cvdata[[i]]$y_predict) |> cbind(mdl)
else
mauc <- rbind(mauc, {MAUC(cvdata[[i]]$y_actual, cvdata[[i]]$y_predict) |> cbind(mdl)})
}
temp <- mauc[, c("mdl", "threshold", "metric")] %>%
pivot_wider(names_from=mdl, values_from=metric, id_cols=threshold)
clr <- c("salmon2", cores[6], "orange2")
mns <- c(
lg=mean(temp$LogReg, na.rm=T),
nb=mean(temp$NBayes, na.rm=T),
rf=mean(temp$RndFst, na.rm=T)) |> round(5)
p1 <- ggplot(data=temp) +
geom_line(aes(x=threshold, y=LogReg), color=clr[1]) +
geom_text(aes(x=.89, y=-.6), label=mns["lg"], color=clr[1]) +
geom_hline(yintercept=0, color=cores[3], linetype="dotted") +
my_ggtheme()
p2 <- ggplot(data=temp) +
geom_line(aes(x=threshold, y=NBayes), color=clr[2]) +
geom_text(aes(x=.88, y=-.6), label=mns["nb"], color=clr[2]) +
geom_hline(yintercept=0, color=cores[3], linetype="dotted") +
my_ggtheme()
p3 <- ggplot(data=temp) +
geom_line(aes(x=threshold, y=RndFst), color=clr[3]) +
geom_text(aes(x=.9, y=.01), label=mns["rf"], color=clr[3]) +
geom_hline(yintercept=0, color=cores[3], linetype="dotted") +
my_ggtheme()
subplot(p1, p2, p3, nrows=3, titleY=T)O melhor resultado que obtivemos foi no Random Forest, já que além de se manter mais alto acima de zero, também o está em mais pontos diferentes. Considerando o espaço entre a curva e a linha de zero como sua área, o Naive Bayes possui uma área maior abaixo da linha que acima, e obtemos um valor negativo, este é um desempenho muito indesejável.
Agora que conseguimos escolher o modelo, vamos definir o limiar de probabilidade que separa o que será classificado como “Acionamento do sinistro” e o que será “Não acionamento do sinistro”. Para isso vamos usar a mesma métrica customizada \(M\), mas desta vez vamos considerar o ponto de limiar que maximiza o valor da métrica. Ainda a partir dos dados da validação cruzada:
temp <- cvdata$RndFst
select_threshold(
temp$y_actual,
temp$y_predict,
plot_=TRUE,
model_name="o Random Forest")Com isto, temos um modelo escolhido e um limiar definido, agora podemos tirar proveito do método holdout para obter métricas de previsão com dados que nunca foram usados para treinar o modelo. Assim poderemos simular uma situação em que estamos colocando-o em produção, exposto à novos dados em tempo real.
Agora chegou a hora de descobrir como este modelo desempenha “no mundo real”, vamos usar o recorte de dados que criamos anteriormente para avaliação para obter uma matriz de confusão:
y <- adf$sinistro
y_pred <- predict(mdl_rf, adf, type="prob")[, "TRUE"]
confusion_matrix({y_pred >= .86}, y)## tp fp tn fn
## 2653 7125 119381 292Baseado nos nossos testes, este é o cenário em que o modelo tem o melhor balanço entre os Positivos Verdadeiros (tp) e os Falsos Negativos (fn). Podemos encontrar o melhor limiar para estes dados e descobrir se o novo valor se diferencia muito do que já obtivemos antes.
select_threshold(y, y_pred, plot_=TRUE, model_name="o Random Forest em produção")Graças à ajuda da validação cruzada conseguimos um valor de limiar muito próximo do ideal para este conjunto de avaliação. Vamos preparar nossa bateria de testes:
f2_score <- function(y_true, y_pred)
FBeta_Score(y_true=y_true, y_pred=y_pred, beta=2)
y_test <- cvdata$RndFst$y_actual
y_test_pred <- cvdata$RndFst$y_predict
tests <- c(`Acurácia`=Accuracy, `Precisão`=Precision, TPR=Recall, F2=f2_score)Para comparar seu desempenho durante o teste, e agora, durante a avaliação:
for (i in seq_along(tests)) {
name <- names(tests)[[i]]
test <- tests[[i]]
tr_result <- test(y_true=y_test, y_pred={y_test_pred >= .86})
ev_result <- test(y_true=y, y_pred={y_pred >= .86})
if (i == 1)
out <- data.frame(
Teste=name,
Treinamento=tr_result,
Avaliação=ev_result
)
else
out <- rbind(out, data.frame(
Teste=name,
Treinamento=tr_result,
Avaliação=ev_result))
}
out_table(out)Teste | Treinamento | Avaliação |
|---|---|---|
Acurácia | 0.9460212 | 0.9427042 |
Precisão | 0.9987133 | 0.9975600 |
TPR | 0.9466018 | 0.9436786 |
F2 | 0.9565844 | 0.9539841 |
Aparentemente não existe muita diferença entre os resultados obtidos no teste e na avaliação, mas percebemos que durante a avaliação, o Random Forest tem desempenho ligeiramente pior, o que está dentro do esperado em qualquer processo de treinamento e avaliação de modelo, e não parece ser uma diferença grande o suficiente para se considerar a possibilidade de overfit.
Obtivemos valores bem altos em todas as métricas, mas isso se deve em boa parte pela “Facilidade” de obter Verdadeiros Negativos (TN), uma vez que o nosso conjunto de dados é desbalanceado. De maneira geral ainda podemos dizer que nosso modelo tem um bom resultado, considerando todas as limitações que enfrentamos.
Obrigado pela leitura!